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A holonomic constraint is used to enforce a constant instantaneous configurational temperature 
on an equilibrium system. Three sets of equations of motion are obtained, differing according to the 
way in which the holonomic constraint is introduced and the phase space distribution function that 
is preserved by the dynamics. Firstly, Gauss' principle of least constraint is used, and it is shown 
that it does not preserve the canonical distribution. Secondly, a modified Haniiltonian is used to 
find a dynamics that provides a restricted microcanonical distribution. Lastly, we provide equations 
that are specifically designed to both satisfy the temperature constraint and produce a restricted 
canonical ensemble. 
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When simulating a molecular system at equilibrium 
one would often like to have it at a constant temperature 
rather than the energy. In a non-equilibrium simulation 
controlling the temperature becomes mandatory as the 
driving forces would otherwise heat up the sytems indef- 
initely. In order to thermostat a system in bulk, away 
from the walls, some synthetic fields can be added to 
the equations of motion. Two important questions arise. 
Firstly, what microscipic measure should one choose for 
the temperature to be controlled? Secondly, what prob- 
ability distribution will the resulting equations generate 
in the system phase space. 

One popular thermostat widely used in molecular dy- 
namics simulations is the isokinetic Gaussian equations 
of motion, which fix the kinetic energy to a desired value 
at all times and generate the isokinetic canonical distri- 
bution Unfortunately, although the kinetic energy 
of molecules is a good measure of the temperature, the 
thermal velocities cannot always be distinguished from 
the streaming velocity which is usually unknown before- 
hand. If the streaming velocity profile is presumed, the 
thermostat will maintain it, acting as a "profile-biased" 
thermostat (PBT), which is often undesirable 0]. 

There are configurational measures for the tempera- 
ture which do not require the peculiar, thermal velocities 
to be known. Several expressions, all equal to one an- 
other in the thermodynamics limit, are known the 
most popular one being (in various notations) 
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where TconF is the configurational temperature in ener- 
getic units, U is the interaction potential energy, g is a 
vector containing the position of all the particles, d is the 
Cartesian dimension of the system, N is the number of 



particles, i labels the degrees of freedom and the angular 
brackets represent an ensemble average. This expression 
has been known as the hypervirial theorem [J and has 
been proved for canonical and microcanonical ensembles 

The first successful implementation of a thermostat 
controlling the configuational temperature was with the 
Nose-Hoover method. The spurious string phase ob- 
served in shear flow simualtions with PBT was eliminated 
when using the configuational thermostat Moreover, 
in the most recent revision of the method Q the pro- 
jected phase space distribution was made canonical. 

In the Nose-Hoover method the dynamics is extended 
by a degree of freedom that describes coupling to the 
thermostat, and the fluctuations of the temperature are 
governed by an extra parameter. It is of interest to see 
if a practical, constant configurational temperature ther- 
mostat can be developed that constrains the original sys- 
tem to a constant instantaneous configurational temper- 
ature. We take Eq.([T]) without the averaging brackets 
as a measure of instantaneous configurational temper- 
ature, by analogy to the Gaussian isokinetic thermostat 
which constrains the instantaneous kinetic energy. In this 
Letter we consider different ways of introducing the con- 
straint into the original equilibrium Hamiltonian equa- 
tions and the effect on the ensemble probability distribu- 
tion. In fact, the resulting equation will be valid not only 
for the configurational thermostat but for any holonomic 
constraint. 

A constraint equation is insufficient to determine the 
constraint reaction forces: an additional principle has to 
be employed in order to find the constrained equations of 
motion [8[ . One can use principles restricting the reaction 
forces, for instance, d'Alambert's principle of orthogonal 
forces or Gauss' principle of least constraint. Alterna- 
tively, there is also the Hamilton variational principle of 
least action. For holonomic constraints, they all lead to 
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the same equations of motion: 



q = -VU - X{q, q)n = -(V;7)|| -n[q 
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where g is a position vector of all the particles, \{q, q) ~ 
q ■ \Vh\~^Vg — n ■ VU is a Lagrange multiplier deter- 
mined as to satisfy the equation of constraint h{q) — ho; 
g{q, q) = h = V/i • q is the conjugate non-holonomic con- 
straint; n — is the unit normal vector to the con- 
straint hypersurface il,h, and (VC/)|| is the component of 
VJ7 tangent to fi/j. In the following we assume unit mass, 
for simplicity. 

One can move to the Hamiltonian formalism by intro- 
ducing the generalized momenta simply as p ^ q. This, 
together with ([2]) give the standard Gaussian equations of 
motion with an isoconfigurational thermostat. The same 
equations are obtained through Dirac's general treatment 
of constrained Hamiltonian systems |9| when the con- 
straint is holonomic: 



q = p 



P = -VC/ - \{q,p)n = (VC/)|| -n\p 
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Let us consider the phase space properties of this dynam- 
ical system. Besides the constrained quantity itself, the 
energy of the system E{q,p) = ^p-p+U is also a constant 
along the actual trajectory: 
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p = 0- A(g,p)(p-n) -0 (4) 



Moreover, if the constraint depends only on intramolec- 
ular distances, as the configurational temperature does, 
then the total momentum P is preserved as well. We 
now investigate the evolution of the phase space elemen- 
tary volumes by calculating the phase space compression 
factor: 

A=— ■ — --0-71— --2n ^ - 
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If we write the Liouville equation with the phase space 
compression factor [l| , 
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the probability distribution density function / is easily 
found and is written over the whole phase space as fol- 
lows, 

/ = V/i • VhS{h - ho)S{g)6{E - Eq)5{P - Pq) (7) 

For some constraints the weight V/i • V/i amounts to a 
constant. This is the case for the bond length constraints, 



for instance. Generally, however, the distribution ^ is 
non-uniform over equal phase space volumes. However, 
the volumes themselves are constantly being stretched 
and compressed by the dynamics as the phase space com- 
pression factor is non-zero (O (albeit being zero on av- 
erage) . If one uses the so called invariant measure of the 
phase space [l^l, which for the dynamics ([3])turns out to 
be just the Vh ■ Vhdqdp [lH, then the density function 
becomes uniform over non-equal but invariant volumes. 
Whichever the interpretation of this weight, the dynam- 
ics will only preserve a non-canonical distribution unless 
V/i ■ V/i is constant, and the values accumulated dur- 
ing the simulation would have to be "reweighted" before 
any comparison with experimental data is possible. It 
is preferable to find such isoconfigurational thermostats 
that would produce a microcanonical, or better still, a 
portion of the Gibbs canonical distribution ~ a restricted 
canonical distribution. 

In order for the configuration to satisfy a holonomic 
constraint, the velocities should be tangent to the con- 
straint hypersurface: q ■ n = 0. Let us take an arbi- 
trary vector p that we can call a peculiar momentum 
and subtract from it its component, which we name a 
"convection" term, perpendicular to the hypersurface: 
q = p — n{p ■ n). We now use the internal energy ex- 
pressed in the coordinates {q,p), 

H = ti + u = P-P-^P-^)\ u (8) 

as a Hamiltonian to generate the following equations of 
motion: 



q = p — n{p ■ n) 

p = -VU + {p ■ n)V{p ■ n) 



(9) 



The Hamiltonian ([5]) is known to be Hamiltonian in re- 
dundant coordinates [l^ . As for any other Hamiltonian 
system, the phase space compression factor is zero, and 
the distribution is the (restricted) microcanonical. 



= 6{h - ho)S{g)S(H - Ho)S{P - Po). 
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Finally, we consider a more general family of equations 
satisfying the holonomic constraint: 



q = p — n{p ■ n) 
p=-(VC/)||+i?(g) 
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where we excluded from the unknown reaction force R a 
term that is to cancel the normal component of the inter- 
molecular interaction force and assumed that R should 
only depend on the configurational information and not 
on the momenta. We will demand the restricted canoni- 
cal distribution. 
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where (3 = I/Tq. According to the LiouviUe equation 
we must then have 
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As we calculate the A and H similarly to ^ and ^ for 
the system ([TT|) . we find the following relation for R: 

—V • n{p ■ n) ~ f3p ■ R 

where the left-hand side can be further modified to see 
the scalar product of p, 

—V • n{p • n) = p • (V • Ti)n 

Since it must be valid for any p we obtain the solution, 

R = -To(V • n)n = -ro[n(V ■ n) + {n ■ V)n] (13) 

where the first term in the brackets is normal and the 
second is tangential to the flh (because n ■ (n ■ y)n — 
i(n- V)(n-n) = 0). 

When TconF = Tq is chosen as the holonomic con- 
straint h{q) = ho, the equations pTjl and p3p describe a 
dynamics that preserves the configurational temperature 
at each time step and generates a canonical distribution 
(fT2)) on the constraint hypersurface. 

In conclusion, three sets of dynamic equation were de- 
rived describing three different isoconfigurational ensem- 
bles: non-canonical jS]), microcanonical (fTO|) : 
and canonical pT|) . p^ . The canonical isoconfigura- 
tional ensemble is deemed to be the best candidate to 
simulate a system at a given temperature that is estab- 
lished instantaneously throughout the system based only 
on the configurational information, which should be es- 
pecially useful in simulation of fiow when the stream- 
ing velocity profile is unknown. Work is underway by 



the authors to numerically test the equations obtained 
and specifically resolve the implementation difficulties, 
such the algebraic unwieldiness and numerical stiffness 
of higher derivatives of potential energy and the choice 
of the starting configuration that corresponds to the de- 
sired temperature. 
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